1 Data

num.splines <- 42
first.columns <- c(
    "speaker",
    "seconds",
    "rec.date",
    "prompt",
    "label",
    "TT.displacement.sm",
    "TT.velocity",
    "TT.velocity.abs",
    "TD.displacement.sm",
    "TD.velocity",
    "TD.velocity.abs"
    )
columns <- c(first.columns,
             paste0(rep(c("X", "Y"), num.splines),
                    "_",
                    rep(1:num.splines, each = 2)
                    )
             )

splines.raw <- list.files(path = "./voicing-effect/results",
                   pattern = "*-aaa.csv",
                   full.names = TRUE) %>%
    map_df(~read.delim(., col.names = columns, na.strings = "*"))
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
splines.closure.raw <- list.files(path = "./voicing-effect/results",
                   pattern = "*-closure.csv",
                   full.names = TRUE) %>%
    map_df(~read.delim(., col.names = columns, na.strings = "*"))
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
splines.phantom.raw <- list.files(path = "./voicing-effect/results",
                   pattern = "*-phantom.csv",
                   full.names = TRUE) %>%
    map_df(~read_tsv(., col_names = columns, na = "*"))
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
rm(num.splines, first.columns, columns)

languages <- read_csv("./voicing-effect/stimuli/languages.csv")
## Parsed with column specification:
## cols(
##   speaker = col_character(),
##   language = col_character()
## )
words <- read_csv("./voicing-effect/stimuli/nonce.csv")
## Parsed with column specification:
## cols(
##   item = col_integer(),
##   word = col_character(),
##   ipa = col_character(),
##   c1 = col_character(),
##   c1phonation = col_character(),
##   vowel = col_character(),
##   anteropost = col_character(),
##   height = col_character(),
##   c2 = col_character(),
##   c2phonation = col_character(),
##   c2place = col_character(),
##   language = col_character()
## )
splines <- splines.raw %>%
    gather(spline, coordinate, matches("[XY]_")) %>%
    separate(spline, c("axis", "fan"), convert = TRUE) %>%
    spread(axis, coordinate) %>%
    mutate(word = word(prompt, 2)) %>%
    left_join(y = languages) %>%
    left_join(y = words) %>%
    mutate_if(is.character, as.factor)
## Joining, by = "speaker"
## Joining, by = c("word", "language")
splines.it <- filter(splines, language == "italian")
splines.pl <- filter(splines, language == "polish")

splines.closure <- splines.closure.raw %>%
    gather(spline, coordinate, matches("[XY]_")) %>%
    separate(spline, c("axis", "fan"), convert = TRUE) %>%
    spread(axis, coordinate) %>%
    mutate(word = word(prompt, 2)) %>%
    left_join(y = languages) %>%
    left_join(y = words) %>%
    mutate_if(is.character, as.factor)
## Joining, by = "speaker"
## Joining, by = c("word", "language")
splines.closure.it <- filter(splines.closure, language == "italian")
splines.closure.pl <- filter(splines.closure, language == "polish")

2 Italian

2.1 Plot at maximum constriction

filter(splines.it, label == "max_TT") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel) +
    scale_colour_discrete(name = "Phonation of C2")
## Warning: Removed 754 rows containing non-finite values (stat_smooth).

filter(splines.it, label == "max_TD") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel) +
    scale_colour_discrete(name = "Phonation of C2")
## Warning: Removed 1127 rows containing non-finite values (stat_smooth).

filter(splines.it, label == "max_TT") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 754 rows containing non-finite values (stat_smooth).

filter(splines.it, label == "max_TD") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 1127 rows containing non-finite values (stat_smooth).

2.2 Plot at peak 1

filter(splines.it, label == "peak1_TT") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 844 rows containing non-finite values (stat_smooth).

filter(splines.it, label == "peak1_TD") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 997 rows containing non-finite values (stat_smooth).

filter(splines.it, label == "peak1_TT") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 844 rows containing non-finite values (stat_smooth).

filter(splines.it, label == "peak1_TD") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 997 rows containing non-finite values (stat_smooth).

2.3 Plot at closure

filter(splines.closure.it, c2place == "coronal") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 937 rows containing non-finite values (stat_smooth).

filter(splines.closure.it, c2place == "velar") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 1150 rows containing non-finite values (stat_smooth).

filter(splines.closure.it, c2place == "coronal") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 937 rows containing non-finite values (stat_smooth).

filter(splines.closure.it, c2place == "velar") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 1150 rows containing non-finite values (stat_smooth).

2.4 GAMMs at maximum constriction

splines.it$c2phonation.ord <- as.ordered(splines.it$c2phonation)
splines.it$c2place.ord <- as.ordered(splines.it$c2place)
splines.it$vowel.ord <- as.ordered(splines.it$vowel)
contrasts(splines.it$c2phonation.ord) <- "contr.treatment"
contrasts(splines.it$c2place.ord) <- "contr.treatment"
contrasts(splines.it$vowel.ord) <- "contr.treatment"
max.it <- filter(splines.it, label %in% c("max_TT", "max_TD"), c1phonation == "voiceless", vowel != "u") %>%
    na.omit() %>%
    arrange(rec.date, fan) %>%
    start_event(column = "fan", event = c("rec.date"))
max.it.gam <- bam(
    Y ~
        c2phonation.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr"),
    data = max.it,
    method = "ML"
)

summary(max.it.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr")
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -5.3889     0.1275 -42.259  < 2e-16 ***
## c2phonation.ord.L   1.2309     0.1809   6.804 1.25e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df       F  p-value    
## s(X)                          8.636  8.938 923.974  < 2e-16 ***
## s(X):c2phonation.ordvoiceless 2.817  3.523   9.283 2.33e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.849   Deviance explained =   85%
## -ML = 8962.5  Scale est. = 43.808    n = 2702
max.it.gam.null <- bam(
    Y ~
        s(X, bs = "cr"),
    data = max.it,
    method = "ML"
)

compareML(max.it.gam, max.it.gam.null)
## max.it.gam: Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr")
## 
## max.it.gam.null: Y ~ s(X, bs = "cr")
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf  Chisq    Df   p.value Sig.
## 1 max.it.gam.null 8999.026   3                            
## 2      max.it.gam 8962.513   6 36.514 3.000 9.592e-16  ***
## 
## AIC difference: -70.19, model max.it.gam has lower AIC.
plot_smooth(max.it.gam, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * X : numeric predictor; with 30 values ranging from -42.260800 to 59.923300.

plot_diff(max.it.gam, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
##  * X : numeric predictor; with 100 values ranging from -42.260800 to 59.923300.

## 
## X window(s) of significant difference(s):
##  -42.260800 - -5.102945
acf_plot(resid(max.it.gam), split_by=list(max.it$rec.date))

max.it.gam.2 <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = max.it,
    method = "fREML"
)

summary(max.it.gam.2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -5.46938    0.08895 -61.490   <2e-16 ***
## c2phonation.ord.L  1.30798    0.12111  10.800   <2e-16 ***
## c2place.ord.L      1.13703    0.12529   9.075   <2e-16 ***
## vowel.ord.L        1.30514    0.12075  10.809   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df       F p-value    
## s(X)                          8.754  8.945 1013.76  <2e-16 ***
## s(X):c2phonation.ordvoiceless 6.292  7.238   18.51  <2e-16 ***
## s(X):c2place.ordvelar         8.701  8.946  314.19  <2e-16 ***
## s(X):vowel.ordo               6.661  7.835   28.53  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.933   Deviance explained = 93.4%
## fREML =   7909  Scale est. = 19.508    n = 2702
max.it.gam.null <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr"),
    data = max.it,
    method = "fREML"
)

compareML(max.it.gam.2, max.it.gam.null)
## max.it.gam.2: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord)
## 
## max.it.gam.null: Y ~ c2phonation.ord + c2place.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr") + s(X, by = c2place.ord, bs = "cr")
## 
## Chi-square test of fREML scores
## -----
##             Model    Score Edf   Chisq    Df  p.value Sig.
## 1 max.it.gam.null 8059.126   9                            
## 2    max.it.gam.2 7908.981  12 150.145 3.000  < 2e-16  ***
## 
## AIC difference: -315.03, model max.it.gam.2 has lower AIC.
plot_smooth(max.it.gam.2, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -42.260800 to 59.923300.

plot_diff(max.it.gam.2, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -42.260800 to 59.923300.

## 
## X window(s) of significant difference(s):
##  -38.132149 - -12.328084
acf_plot(resid(max.it.gam.2), split_by=list(max.it$rec.date))

max.it.gamm <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr"),
    data = max.it,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(max.it.gamm)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr", 
##     m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X, 
##     by = vowel.ord, bs = "cr")
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)   
## (Intercept)        -9.8210     4.4120  -2.226  0.02611 * 
## c2phonation.ord.L   1.4489     0.4473   3.239  0.00121 **
## c2place.ord.L       0.8747     0.4548   1.923  0.05454 . 
## vowel.ord.L         1.2193     0.4522   2.696  0.00706 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                   edf  Ref.df       F  p-value    
## s(X)                            6.161   6.763   6.500 1.21e-07 ***
## s(X):c2phonation.ordvoiceless   5.441   6.330   4.058 0.000408 ***
## s(X,rec.date)                 234.516 412.000   4.536  < 2e-16 ***
## s(X,speaker)                    6.481   8.000  32.197  < 2e-16 ***
## s(X):c2place.ordvelar           8.811   8.953 197.900  < 2e-16 ***
## s(X):vowel.ordo                 6.656   7.497  11.144 1.54e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.965   Deviance explained = 96.9%
## fREML = 7334.8  Scale est. = 10.057    n = 2702
max.it.gam.null <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = max.it,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(max.it.gamm, max.it.gam.null)
## max.it.gamm: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr", 
##     m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X, 
##     by = vowel.ord, bs = "cr")
## 
## max.it.gam.null: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord)
## 
## Chi-square test of fREML scores
## -----
##             Model    Score Edf  Chisq    Df  p.value Sig.
## 1 max.it.gam.null 7430.089  15                           
## 2     max.it.gamm 7334.843  18 95.246 3.000  < 2e-16  ***
## 
## AIC difference: -141.27, model max.it.gamm has lower AIC.
plot_smooth(max.it.gamm, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -42.260800 to 59.923300. 
##  * rec.date : factor; set to the value(s): 29/11/2016 15:11:03. 
##  * speaker : factor; set to the value(s): sc01.

plot_diff(max.it.gamm, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -42.260800 to 59.923300. 
##  * rec.date : factor; set to the value(s): 29/11/2016 15:11:03. 
##  * speaker : factor; set to the value(s): sc01.

## 
## X window(s) of significant difference(s):
##  -42.260800 - -15.424572
acf_plot(resid(max.it.gamm), split_by=list(max.it$rec.date))

2.4.1 AR1 model

r1 <- start_value_rho(max.it.gamm)

max.it.gam.AR <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr"),
    data = max.it,
    method = "fREML",
    rho = r1,
    AR.start = max.it$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(max.it.gam.AR)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr", 
##     m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X, 
##     by = vowel.ord, bs = "cr")
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -7.4983     3.6253  -2.068  0.03871 *  
## c2phonation.ord.L   1.4251     0.3367   4.233 2.39e-05 ***
## c2place.ord.L       1.0267     0.3385   3.033  0.00244 ** 
## vowel.ord.L         1.0714     0.3364   3.184  0.00147 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                   edf  Ref.df       F  p-value    
## s(X)                            7.254   7.812   7.852 4.16e-09 ***
## s(X):c2phonation.ordvoiceless   7.172   8.020   7.913 1.42e-10 ***
## s(X,rec.date)                 134.559 412.000   0.807  < 2e-16 ***
## s(X,speaker)                    5.706   8.000  25.634  < 2e-16 ***
## s(X):c2place.ordvelar           8.709   8.936 135.867  < 2e-16 ***
## s(X):vowel.ordo                 5.644   6.654   8.406 1.41e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.953   Deviance explained = 95.6%
## fREML = 5823.8  Scale est. = 7.9222    n = 2702
max.it.gam.AR.null <- bam(
    Y ~
#        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr"),
    data = max.it,
    method = "fREML",
    rho = r1,
    AR.start = max.it$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(max.it.gam.AR.null, max.it.gam.AR)
## max.it.gam.AR.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, rec.date, 
##     bs = "fs", xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr")
## 
## max.it.gam.AR: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr", 
##     m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X, 
##     by = vowel.ord, bs = "cr")
## 
## Chi-square test of fREML scores
## -----
##                Model    Score Edf  Chisq    Df   p.value Sig.
## 1 max.it.gam.AR.null 5861.708  15                            
## 2      max.it.gam.AR 5823.848  18 37.860 3.000 2.539e-16  ***
## Warning in compareML(max.it.gam.AR.null, max.it.gam.AR): AIC is not
## reliable, because an AR1 model is included (rho1 = 0.725285, rho2 =
## 0.725285).
acf_plot(resid(max.it.gam.AR), split_by=list(max.it$rec.date))

acf_resid(max.it.gam.AR, split_pred = "AR.start")

plot_smooth(max.it.gam.AR, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -42.260800 to 59.923300. 
##  * rec.date : factor; set to the value(s): 29/11/2016 15:11:03. 
##  * speaker : factor; set to the value(s): sc01.

plot_diff(max.it.gam.AR, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -42.260800 to 59.923300. 
##  * rec.date : factor; set to the value(s): 29/11/2016 15:11:03. 
##  * speaker : factor; set to the value(s): sc01.

## 
## X window(s) of significant difference(s):
##  -41.228637 - -17.488897

2.5 GAMMs at closure

splines.closure.it$c2phonation.ord <- as.ordered(splines.closure.it$c2phonation)
splines.closure.it$c2place.ord <- as.ordered(splines.closure.it$c2place)
splines.closure.it$vowel.ord <- as.ordered(splines.closure.it$vowel)
contrasts(splines.closure.it$c2phonation.ord) <- "contr.treatment"
contrasts(splines.closure.it$c2place.ord) <- "contr.treatment"
contrasts(splines.closure.it$vowel.ord) <- "contr.treatment"
closure.it <- splines.closure.it %>%
    filter(c1phonation == "voiceless", vowel != "u") %>%
    na.omit() %>%
    arrange(rec.date, fan) %>%
    start_event(column = "fan", event = c("rec.date"))
closure.it.gam <- bam(
    Y ~
        c2phonation.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr"),
    data = closure.it,
    method = "ML"
)

summary(closure.it.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr")
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -5.9307     0.1124 -52.776  < 2e-16 ***
## c2phonation.ord.L   0.7339     0.1591   4.614 4.14e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df        F  p-value    
## s(X)                          8.676  8.927 1037.842  < 2e-16 ***
## s(X):c2phonation.ordvoiceless 6.166  7.122    9.889 2.09e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.86   Deviance explained = 86.1%
## -ML = 9003.8  Scale est. = 35.076    n = 2806
closure.it.gam.null <- bam(
    Y ~
        s(X, bs = "cr"),
    data = closure.it,
    method = "ML"
)

compareML(closure.it.gam, closure.it.gam.null)
## closure.it.gam: Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr")
## 
## closure.it.gam.null: Y ~ s(X, bs = "cr")
## 
## Chi-square test of ML scores
## -----
##                 Model    Score Edf  Chisq    Df   p.value Sig.
## 1 closure.it.gam.null 9039.474   3                            
## 2      closure.it.gam 9003.765   6 35.709 3.000 2.120e-15  ***
## 
## AIC difference: -81.83, model closure.it.gam has lower AIC.
plot_smooth(closure.it.gam, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * X : numeric predictor; with 30 values ranging from -45.980000 to 64.280100.

plot_diff(closure.it.gam, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
##  * X : numeric predictor; with 100 values ranging from -45.980000 to 64.280100.

## 
## X window(s) of significant difference(s):
##  -35.956355 - -19.250279
##  -5.885418 - 10.820658
##  27.526733 - 38.664117
acf_plot(resid(closure.it.gam), split_by=list(closure.it$rec.date))

closure.it.gam.2 <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr"),
    data = closure.it,
    method = "fREML"
)

summary(closure.it.gam.2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord, bs = "cr")
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -5.90863    0.08051 -73.390  < 2e-16 ***
## c2phonation.ord.L  0.94698    0.10990   8.617  < 2e-16 ***
## c2place.ord.L      0.73290    0.11296   6.488 1.03e-10 ***
## vowel.ord.L        1.29048    0.11020  11.710  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df       F p-value    
## s(X)                          8.749  8.927 1109.85  <2e-16 ***
## s(X):c2phonation.ordvoiceless 7.178  8.029   23.34  <2e-16 ***
## s(X):c2place.ordvelar         8.736  8.959  277.75  <2e-16 ***
## s(X):vowel.ordo               7.730  8.440   48.76  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.934   Deviance explained = 93.5%
## fREML = 7989.8  Scale est. = 16.566    n = 2806
closure.it.gam.null <- bam(
    Y ~
#        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr"),
    data = closure.it,
    method = "fREML"
)

compareML(closure.it.gam.2, closure.it.gam.null)
## closure.it.gam.2: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord, bs = "cr")
## 
## closure.it.gam.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord, bs = "cr")
## 
## Chi-square test of fREML scores
## -----
##                 Model    Score Edf   Chisq    Df  p.value Sig.
## 1 closure.it.gam.null 8102.335   9                            
## 2    closure.it.gam.2 7989.848  12 112.487 3.000  < 2e-16  ***
## 
## AIC difference: -242.02, model closure.it.gam.2 has lower AIC.
plot_smooth(closure.it.gam.2, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -45.980000 to 64.280100.

plot_diff(closure.it.gam.2, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -45.980000 to 64.280100.

## 
## X window(s) of significant difference(s):
##  -45.980000 - -43.752523
##  -35.956355 - -17.022802
##  -2.544203 - 16.389349
acf_plot(resid(closure.it.gam.2), split_by=list(closure.it$rec.date))

closure.it.gamm <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr"),
    data = closure.it,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(closure.it.gamm)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr", 
##     m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X, 
##     by = vowel.ord, bs = "cr")
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -8.8567     4.2518  -2.083  0.03735 *  
## c2phonation.ord.L   1.0214     0.4111   2.485  0.01303 *  
## c2place.ord.L       0.6041     0.4129   1.463  0.14358    
## vowel.ord.L         1.3767     0.4126   3.337  0.00086 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                   edf  Ref.df       F  p-value    
## s(X)                            6.835   7.419   9.060 1.89e-11 ***
## s(X):c2phonation.ordvoiceless   7.272   8.048   8.385 2.50e-11 ***
## s(X,rec.date)                 249.926 402.000   5.429  < 2e-16 ***
## s(X,speaker)                    5.928   8.000  33.672  < 2e-16 ***
## s(X):c2place.ordvelar           8.801   8.958 205.734  < 2e-16 ***
## s(X):vowel.ordo                 8.593   8.884  23.072  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.973   Deviance explained = 97.5%
## fREML = 7121.2  Scale est. = 6.9052    n = 2806
closure.it.gamm.null <- bam(
    Y ~
#        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr"),
    data = closure.it,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(closure.it.gamm, closure.it.gamm.null)
## closure.it.gamm: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr", 
##     m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X, 
##     by = vowel.ord, bs = "cr")
## 
## closure.it.gamm.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, rec.date, 
##     bs = "fs", xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr")
## 
## Chi-square test of fREML scores
## -----
##                  Model    Score Edf  Chisq    Df   p.value Sig.
## 1 closure.it.gamm.null 7157.603  15                            
## 2      closure.it.gamm 7121.171  18 36.432 3.000 1.040e-15  ***
## 
## AIC difference: 31.42, model closure.it.gamm.null has lower AIC.
AIC(closure.it.gamm, closure.it.gamm.null)
##                            df      AIC
## closure.it.gamm      295.5252 13668.49
## closure.it.gamm.null 311.7325 13637.07
plot_smooth(closure.it.gamm, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -45.980000 to 64.280100. 
##  * rec.date : factor; set to the value(s): 12/12/2016 14:44:52. 
##  * speaker : factor; set to the value(s): sc01.

plot_diff(closure.it.gamm, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -45.980000 to 64.280100. 
##  * rec.date : factor; set to the value(s): 12/12/2016 14:44:52. 
##  * speaker : factor; set to the value(s): sc01.

## 
## X window(s) of significant difference(s):
##  -40.411308 - -18.136540
acf_plot(resid(closure.it.gamm), split_by=list(closure.it$rec.date))

2.5.1 AR1 model

r1 <- start_value_rho(closure.it.gamm)

closure.it.gam.AR <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr"),
    data = closure.it,
    method = "fREML",
    rho = r1,
    AR.start = closure.it$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(closure.it.gam.AR)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr", 
##     m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X, 
##     by = vowel.ord, bs = "cr")
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -8.1777     3.8108  -2.146 0.031970 *  
## c2phonation.ord.L   1.0379     0.3333   3.114 0.001867 ** 
## c2place.ord.L       0.7008     0.3345   2.095 0.036297 *  
## vowel.ord.L         1.2586     0.3341   3.768 0.000168 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                   edf  Ref.df       F  p-value    
## s(X)                            7.185   7.753   7.214 4.45e-09 ***
## s(X):c2phonation.ordvoiceless   6.983   7.901   8.287 5.12e-11 ***
## s(X,rec.date)                 186.323 402.000   1.730  < 2e-16 ***
## s(X,speaker)                    5.740   8.000  32.534  < 2e-16 ***
## s(X):c2place.ordvelar           8.730   8.946 144.770  < 2e-16 ***
## s(X):vowel.ordo                 7.860   8.526  17.379  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.965   Deviance explained = 96.8%
## fREML = 5533.1  Scale est. = 5.2451    n = 2806
closure.it.gam.AR.null <- bam(
    Y ~
#        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr"),
    data = closure.it,
    method = "fREML",
    rho = r1,
    AR.start = closure.it$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(closure.it.gam.AR, closure.it.gam.AR.null)
## closure.it.gam.AR: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr", 
##     m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X, 
##     by = vowel.ord, bs = "cr")
## 
## closure.it.gam.AR.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, rec.date, 
##     bs = "fs", xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr")
## 
## Chi-square test of fREML scores
## -----
##                    Model    Score Edf  Chisq    Df   p.value Sig.
## 1 closure.it.gam.AR.null 5566.155  15                            
## 2      closure.it.gam.AR 5533.075  18 33.079 3.000 2.835e-14  ***
## Warning in compareML(closure.it.gam.AR, closure.it.gam.AR.null): AIC is
## not reliable, because an AR1 model is included (rho1 = 0.725722, rho2 =
## 0.725722).
acf_plot(resid(closure.it.gam.AR), split_by=list(closure.it$rec.date))

acf_resid(closure.it.gam.AR, split_pred = "AR.start")

plot_smooth(closure.it.gam.AR, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -45.980000 to 64.280100. 
##  * rec.date : factor; set to the value(s): 12/12/2016 14:44:52. 
##  * speaker : factor; set to the value(s): sc01.

plot_diff(closure.it.gam.AR, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -45.980000 to 64.280100. 
##  * rec.date : factor; set to the value(s): 12/12/2016 14:44:52. 
##  * speaker : factor; set to the value(s): sc01.

## 
## X window(s) of significant difference(s):
##  -40.411308 - -15.909064

2.6 GAMMs at peak 1

splines.it$c2phonation.ord <- as.ordered(splines.it$c2phonation)
contrasts(splines.it$c2phonation.ord) <- "contr.treatment"
peak.it <- filter(splines.it, label %in% c("peak1_TT", "peak1_TD"))

peak.it.gam <- bam(
    Y ~
        c2phonation.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr"),
    data = peak.it,
    method = "ML"
)

summary(peak.it.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr")
## 
## Parametric coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -5.8024     0.1809 -32.083   <2e-16 ***
## c2phonation.ordvoiceless   0.5019     0.2541   1.976   0.0482 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df       F p-value    
## s(X)                          8.691  8.958 562.990  <2e-16 ***
## s(X):c2phonation.ordvoiceless 1.695  2.137   0.932   0.425    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.59   Deviance explained =   59%
## -ML =  22127  Scale est. = 96.198    n = 5971
peak.it.gam.null <- bam(
    Y ~
        s(X, bs = "cr"),
    data = peak.it,
    method = "ML"
)

compareML(peak.it.gam, peak.it.gam.null)
## peak.it.gam: Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr")
## 
## peak.it.gam.null: Y ~ s(X, bs = "cr")
## 
## Chi-square test of ML scores
## -----
##              Model    Score Edf Chisq    Df p.value Sig.
## 1 peak.it.gam.null 22129.00   3                         
## 2      peak.it.gam 22126.83   6 2.171 3.000   0.227     
## 
## AIC difference: -0.63, model peak.it.gam has lower AIC.
## Warning in compareML(peak.it.gam, peak.it.gam.null): Only small difference in ML...
plot_smooth(peak.it.gam, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * X : numeric predictor; with 30 values ranging from -48.091200 to 64.315100.

plot_diff(peak.it.gam, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
##  * X : numeric predictor; with 100 values ranging from -48.091200 to 64.315100.

## 
## Difference is not significant.

2.7 Comparison of voiced consonants at closure and maximum constriction

voiced.it <- rbind(max.it, closure.it) %>%
    filter(c2phonation == "voiced") %>%
    mutate(label.ord = ifelse(label %in% c("max_TT", "max_TD"), "maximum", "closure"))

voiced.it$label.ord <- as.ordered(voiced.it$label.ord)
contrasts(voiced.it$label.ord) <- "contr.treatment"
filter(voiced.it, c2place == "coronal") %>%
    ggplot(aes(X, Y, colour = label.ord)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)

voiced.it.gamm <- bam(
    Y ~
        label.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = label.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr"),
    data = voiced.it,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(voiced.it.gamm)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ label.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, 
##     by = label.ord, bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", 
##     m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr", m = 1, 
##     k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr")
## 
## Parametric coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -10.3054     4.7990  -2.147   0.0319 *  
## label.ordmaximum  -1.3368     0.1183 -11.303   <2e-16 ***
## c2place.ord.L      0.9681     0.7984   1.212   0.2255    
## vowel.ord.L        1.3332     0.7917   1.684   0.0923 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                           edf  Ref.df      F  p-value    
## s(X)                    5.881   6.501  7.462 1.32e-08 ***
## s(X):label.ordmaximum   7.188   8.018 69.108  < 2e-16 ***
## s(X,rec.date)         154.503 209.000  8.621  < 2e-16 ***
## s(X,speaker)            6.652   8.000 27.102  < 2e-16 ***
## s(X):c2place.ordvelar   8.742   8.929 90.637  < 2e-16 ***
## s(X):vowel.ordo         7.535   8.218 12.681  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.966   Deviance explained = 96.8%
## fREML = 7583.5  Scale est. = 9.5337    n = 2846
r1 <- start_value_rho(voiced.it.gamm)

voiced.it.gam.AR <- bam(
    Y ~
        label.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = label.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr"),
    data = voiced.it,
    method = "fREML",
    rho = r1,
    AR.start = voiced.it$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(voiced.it.gam.AR)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ label.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, 
##     by = label.ord, bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", 
##     m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr", m = 1, 
##     k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr")
## 
## Parametric coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -10.4446     4.2142  -2.478   0.0133 *  
## label.ordmaximum  -1.2193     0.2554  -4.774  1.9e-06 ***
## c2place.ord.L      0.9768     0.4913   1.988   0.0469 *  
## vowel.ord.L        1.0384     0.4902   2.118   0.0342 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                           edf  Ref.df       F  p-value    
## s(X)                    5.901   6.590   6.529 1.76e-07 ***
## s(X):label.ordmaximum   7.228   8.101  17.534  < 2e-16 ***
## s(X,rec.date)         109.376 209.000   2.344  < 2e-16 ***
## s(X,speaker)            6.672   8.000  25.111  < 2e-16 ***
## s(X):c2place.ordvelar   8.673   8.924 125.385  < 2e-16 ***
## s(X):vowel.ordo         7.029   7.916  13.625  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.957   Deviance explained = 95.9%
## fREML = 5914.4  Scale est. = 7.3839    n = 2846
plot_smooth(voiced.it.gam.AR, view = "X",
            plot_all= "label.ord", rug = FALSE)
## Summary:
##  * label.ord : factor; set to the value(s): closure, maximum. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -42.852800 to 64.280100. 
##  * rec.date : factor; set to the value(s): 29/11/2016 15:21:30. 
##  * speaker : factor; set to the value(s): sc01.

plot_diff(voiced.it.gam.AR, view = "X",
          comp = list(label.ord = c("maximum", "closure")))
## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -42.852800 to 64.280100. 
##  * rec.date : factor; set to the value(s): 29/11/2016 15:21:30. 
##  * speaker : factor; set to the value(s): sc01.

## 
## X window(s) of significant difference(s):
##  -38.524198 - -15.799037
##  -2.813231 - 36.144187
acf_resid(voiced.it.gam.AR, split_pred = "AR.start")

3 Polish

3.1 Plot at maximum constriction

filter(splines.pl, label == "max_TT") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1.5) +
    facet_grid(speaker ~ vowel) +
    scale_colour_discrete(name = "Phonation of C2")
## Warning: Removed 274 rows containing non-finite values (stat_smooth).

filter(splines.pl, label == "max_TD") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel) +
    scale_colour_discrete(name = "Phonation of C2")
## Warning: Removed 328 rows containing non-finite values (stat_smooth).

filter(splines.pl, label == "max_TT") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1.5) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 274 rows containing non-finite values (stat_smooth).

filter(splines.pl, label == "max_TD") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1.5) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 328 rows containing non-finite values (stat_smooth).

3.2 Plot at peak 1

filter(splines.pl, label == "peak1_TT") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1.5) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 273 rows containing non-finite values (stat_smooth).

filter(splines.pl, label == "peak1_TD") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1.5) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 302 rows containing non-finite values (stat_smooth).

filter(splines.pl, label == "peak1_TT") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1.5) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 273 rows containing non-finite values (stat_smooth).

filter(splines.pl, label == "peak1_TD") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1.5) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 302 rows containing non-finite values (stat_smooth).

3.3 Plot at closure

filter(splines.closure.pl, c2place == "coronal") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1.5) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 247 rows containing non-finite values (stat_smooth).

filter(splines.closure.pl, c2place == "velar") %>%
    ggplot(aes(X, Y, colour = c2phonation, group = rec.date)) +
    geom_line(stat = "smooth", method = "loess", alpha = 0.5) +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 284 rows containing non-finite values (stat_smooth).

filter(splines.closure.pl, c2place == "coronal") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 247 rows containing non-finite values (stat_smooth).

filter(splines.closure.pl, c2place == "velar") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_line(stat = "smooth", method = "loess") +
    coord_fixed(ratio = 1) +
    facet_grid(speaker ~ vowel)
## Warning: Removed 284 rows containing non-finite values (stat_smooth).

3.4 GAMMs at maximum constriction ab03

splines.pl$c2phonation.ord <- as.ordered(splines.pl$c2phonation)
splines.pl$c2place.ord <- as.ordered(splines.pl$c2place)
splines.pl$vowel.ord <- as.ordered(splines.pl$vowel)
contrasts(splines.pl$c2phonation.ord) <- "contr.treatment"
contrasts(splines.pl$c2place.ord) <- "contr.treatment"
contrasts(splines.pl$vowel.ord) <- "contr.treatment"
max.pl <- splines.pl %>%
    filter(label %in% c("max_TT", "max_TD"),
           speaker == "ab03",
           vowel != "u"
    ) %>%
    na.omit() %>%
    arrange(rec.date, fan) %>%
    start_event(column = "fan", event = c("rec.date"))
max.pl.gam <- bam(
    Y ~
        c2phonation.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr"),
    data = max.pl,
    method = "ML"
)

summary(max.pl.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr")
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.9310     0.1236   7.530 8.09e-14 ***
## c2phonation.ord.L   0.7092     0.1752   4.048 5.40e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df       F p-value    
## s(X)                          8.865  8.993 332.184  <2e-16 ***
## s(X):c2phonation.ordvoiceless 2.129  2.662   3.826  0.0208 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.705   Deviance explained = 70.7%
## -ML = 5426.3  Scale est. = 26.88     n = 1763
max.pl.gam.null <- bam(
    Y ~
        s(X, bs = "cr"),
    data = max.pl,
    method = "ML"
)

compareML(max.pl.gam, max.pl.gam.null)
## max.pl.gam: Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr")
## 
## max.pl.gam.null: Y ~ s(X, bs = "cr")
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf  Chisq    Df   p.value Sig.
## 1 max.pl.gam.null 5437.967   3                            
## 2      max.pl.gam 5426.306   6 11.661 3.000 3.461e-05  ***
## 
## AIC difference: -20.14, model max.pl.gam has lower AIC.
plot_smooth(max.pl.gam, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * X : numeric predictor; with 30 values ranging from -34.964200 to 70.078500.

plot_diff(max.pl.gam, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
##  * X : numeric predictor; with 100 values ranging from -34.964200 to 70.078500.

## 
## X window(s) of significant difference(s):
##  23.392856 - 70.078500
max.pl.gam.2 <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = max.pl,
    method = "fREML"
)

summary(max.pl.gam.2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0.86975    0.08813   9.869  < 2e-16 ***
## c2phonation.ord.L  0.63403    0.12287   5.160 2.75e-07 ***
## c2place.ord.L     -0.17665    0.12468  -1.417    0.157    
## vowel.ord.L       -0.10397    0.12277  -0.847    0.397    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df      F p-value    
## s(X)                          8.660  8.929 204.37  <2e-16 ***
## s(X):c2phonation.ordvoiceless 2.511  3.135   3.37  0.0168 *  
## s(X):c2place.ordvelar         8.412  8.879 174.24  <2e-16 ***
## s(X):vowel.ordo               8.275  8.848  32.24  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.856   Deviance explained = 85.8%
## fREML = 4825.9  Scale est. = 13.168    n = 1763
max.pl.gam.null <- bam(
    Y ~
#        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = max.pl,
    method = "fREML"
)

compareML(max.pl.gam.2, max.pl.gam.null)
## max.pl.gam.2: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord)
## 
## max.pl.gam.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord)
## 
## Chi-square test of fREML scores
## -----
##             Model    Score Edf  Chisq    Df   p.value Sig.
## 1 max.pl.gam.null 4842.954   9                            
## 2    max.pl.gam.2 4825.858  12 17.095 3.000 1.806e-07  ***
## 
## AIC difference: -31.01, model max.pl.gam.2 has lower AIC.
plot_smooth(max.pl.gam.2, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -34.964200 to 70.078500.

plot_diff(max.pl.gam.2, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -34.964200 to 70.078500.

## 
## X window(s) of significant difference(s):
##  -30.720051 - -11.621378
##  20.209743 - 70.078500
max.pl.gamm <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
#        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr"),
    data = max.pl,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(max.pl.gamm)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr")
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)
## (Intercept)         0.8847     0.5885   1.503    0.133
## c2phonation.ord.L   1.2735     0.7781   1.637    0.102
## c2place.ord.L      -0.2137     0.8299  -0.257    0.797
## vowel.ord.L        -0.1053     0.8298  -0.127    0.899
## 
## Approximate significance of smooth terms:
##                                   edf  Ref.df      F p-value    
## s(X)                            8.536   8.832 49.903 < 2e-16 ***
## s(X):c2phonation.ordvoiceless   6.022   6.993  2.801 0.00728 ** 
## s(X,rec.date)                 202.133 227.000 26.452 < 2e-16 ***
## s(X):c2place.ordvelar           8.749   8.933 30.231 < 2e-16 ***
## s(X):vowel.ordo                 8.748   8.931 18.211 < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.968   Deviance explained = 97.3%
## fREML = 3926.3  Scale est. = 2.8886    n = 1763
max.pl.gamm.null <- bam(
    Y ~
#        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
#        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr"),
    data = max.pl,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(max.pl.gamm, max.pl.gam.null)
## max.pl.gamm: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr")
## 
## max.pl.gam.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord)
## 
## Chi-square test of fREML scores
## -----
##             Model    Score Edf   Chisq    Df  p.value Sig.
## 1 max.pl.gam.null 4842.954   9                            
## 2     max.pl.gamm 3926.300  15 916.653 6.000  < 2e-16  ***
## 
## AIC difference: -2517.12, model max.pl.gamm has lower AIC.
plot_smooth(max.pl.gamm, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -34.964200 to 70.078500. 
##  * rec.date : factor; set to the value(s): 24/03/2017 14:44:40.

plot_diff(max.pl.gamm, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -34.964200 to 70.078500. 
##  * rec.date : factor; set to the value(s): 24/03/2017 14:44:40.

## 
## X window(s) of significant difference(s):
##  -34.964200 - -22.231752

3.4.1 AR1 model

r1 <- start_value_rho(max.pl.gamm)

max.pl.gam.AR <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr"),
    data = max.pl,
    method = "fREML",
    rho = r1,
    AR.start = max.pl$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(max.pl.gam.AR)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr")
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)  
## (Intercept)         0.8682     0.5156   1.684   0.0924 .
## c2phonation.ord.L   1.1093     0.6731   1.648   0.0995 .
## c2place.ord.L      -0.1025     0.7273  -0.141   0.8879  
## vowel.ord.L        -0.1336     0.7276  -0.184   0.8543  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                   edf  Ref.df      F p-value    
## s(X)                            8.456   8.824 43.295  <2e-16 ***
## s(X):c2phonation.ordvoiceless   4.942   5.993  1.600   0.124    
## s(X,rec.date)                 187.491 227.000  8.812  <2e-16 ***
## s(X):c2place.ordvelar           8.672   8.920 27.668  <2e-16 ***
## s(X):vowel.ordo                 8.729   8.936 15.496  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.962   Deviance explained = 96.7%
## fREML = 3256.1  Scale est. = 2.4815    n = 1763
max.pl.gam.AR.null <- bam(
    Y ~
#        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr"),
    data = max.pl,
    method = "fREML",
    rho = r1,
    AR.start = max.pl$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(max.pl.gam.AR.null, max.pl.gam.AR)
## max.pl.gam.AR.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, rec.date, 
##     bs = "fs", xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord, bs = "cr")
## 
## max.pl.gam.AR: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr")
## 
## Chi-square test of fREML scores
## -----
##                Model    Score Edf Chisq    Df   p.value Sig.
## 1 max.pl.gam.AR.null 3265.122  12                           
## 2      max.pl.gam.AR 3256.061  15 9.062 3.000 4.148e-04  ***
## Warning in compareML(max.pl.gam.AR.null, max.pl.gam.AR): AIC is not
## reliable, because an AR1 model is included (rho1 = 0.615458, rho2 =
## 0.615458).
acf_plot(resid(max.pl.gam.AR), split_by=list(max.pl$rec.date))

acf_resid(max.pl.gam.AR, split_pred = "AR.start")

plot_smooth(max.pl.gam.AR, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -34.964200 to 70.078500. 
##  * rec.date : factor; set to the value(s): 24/03/2017 14:44:40.

plot_diff(max.pl.gam.AR, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -34.964200 to 70.078500. 
##  * rec.date : factor; set to the value(s): 24/03/2017 14:44:40.

## 
## X window(s) of significant difference(s):
##  50.979827 - 54.162939

3.5 GAMMs at maximum constriction ps02

max.pl.ps02 <- splines.pl %>%
    filter(label %in% c("max_TT", "max_TD"),
           speaker == "ps02",
           vowel != "u",
           X < 40,
           X > -15
    ) %>%
    na.omit() %>%
    arrange(rec.date, fan) %>%
    start_event(column = "fan", event = c("rec.date"))
max.pl.gamm.ps02 <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
#        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr"),
    data = max.pl.ps02,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(max.pl.gamm.ps02)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr")
## 
## Parametric coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -13.14685    0.18610 -70.644  < 2e-16 ***
## c2phonation.ord.L  -0.08762    0.26224  -0.334 0.738325    
## c2place.ord.L       1.83074    0.26315   6.957  4.8e-12 ***
## vowel.ord.L         0.96946    0.26082   3.717 0.000208 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                   edf  Ref.df      F p-value    
## s(X)                            8.207   8.665 64.736 < 2e-16 ***
## s(X):c2phonation.ordvoiceless   4.291   5.194  2.243 0.04978 *  
## s(X,rec.date)                 144.962 292.000  2.746 < 2e-16 ***
## s(X):c2place.ordvelar           8.804   8.951 90.910 < 2e-16 ***
## s(X):vowel.ordo                 2.379   2.883  5.123 0.00166 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.863   Deviance explained = 87.5%
## fREML = 5726.5  Scale est. = 13.737    n = 2027

3.5.1 AR1 model

r1 <- start_value_rho(max.pl.gamm.ps02)

max.pl.gam.ps02.AR <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr"),
    data = max.pl.ps02,
    method = "fREML",
    rho = r1,
    AR.start = max.pl.ps02$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(max.pl.gam.ps02.AR)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr")
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -13.6822     0.2054 -66.597  < 2e-16 ***
## c2phonation.ord.L   0.1534     0.2793   0.549 0.582947    
## c2place.ord.L       1.6141     0.2904   5.559  3.1e-08 ***
## vowel.ord.L         0.8489     0.2306   3.681 0.000239 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                   edf  Ref.df      F  p-value    
## s(X)                            8.163   8.667 65.792  < 2e-16 ***
## s(X):c2phonation.ordvoiceless   4.255   5.311  1.418 0.219832    
## s(X,rec.date)                 125.209 292.000  1.364  < 2e-16 ***
## s(X):c2place.ordvelar           8.602   8.883 41.201  < 2e-16 ***
## s(X):vowel.ordo                 1.000   1.001 11.903 0.000572 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.828   Deviance explained = 83.8%
## fREML = 4344.2  Scale est. = 7.8385    n = 2027
max.pl.gam.ps02.AR.null <- bam(
    Y ~
#        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr"),
    data = max.pl.ps02,
    method = "fREML",
    rho = r1,
    AR.start = max.pl.ps02$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(max.pl.gam.ps02.AR.null, max.pl.gam.ps02.AR)
## max.pl.gam.ps02.AR.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, rec.date, 
##     bs = "fs", xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord, bs = "cr")
## 
## max.pl.gam.ps02.AR: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr")
## 
## Chi-square test of fREML scores
## -----
##                     Model    Score Edf Chisq    Df p.value Sig.
## 1 max.pl.gam.ps02.AR.null 4347.654  12                         
## 2      max.pl.gam.ps02.AR 4344.226  15 3.428 3.000   0.077
## Warning in compareML(max.pl.gam.ps02.AR.null, max.pl.gam.ps02.AR): AIC is
## not reliable, because an AR1 model is included (rho1 = 0.733987, rho2 =
## 0.733987).
## Warning in compareML(max.pl.gam.ps02.AR.null, max.pl.gam.ps02.AR): Only small difference in fREML...
acf_plot(resid(max.pl.gam.ps02.AR), split_by=list(max.pl.ps02$rec.date))

acf_resid(max.pl.gam.ps02.AR, split_pred = "AR.start")

plot_smooth(max.pl.gam.ps02.AR, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): velar. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -14.986000 to 39.996900. 
##  * rec.date : factor; set to the value(s): 07/02/2017 16:30:56.

plot_diff(max.pl.gam.ps02.AR, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
##  * c2place.ord : factor; set to the value(s): velar. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -14.986000 to 39.996900. 
##  * rec.date : factor; set to the value(s): 07/02/2017 16:30:56.

## 
## Difference is not significant.

3.6 GAMMs at closure

splines.closure.pl$c2phonation.ord <- as.ordered(splines.closure.pl$c2phonation)
splines.closure.pl$c2place.ord <- as.ordered(splines.closure.pl$c2place)
splines.closure.pl$vowel.ord <- as.ordered(splines.closure.pl$vowel)
contrasts(splines.closure.pl$c2phonation.ord) <- "contr.treatment"
contrasts(splines.closure.pl$c2place.ord) <- "contr.treatment"
contrasts(splines.closure.pl$vowel.ord) <- "contr.treatment"
closure.pl <- splines.closure.pl %>%
    filter(speaker == "ab03", vowel != "u", X > -15) %>%
    na.omit() %>%
    arrange(rec.date, fan) %>%
    start_event(column = "fan", event = c("rec.date"))
closure.pl.gam <- bam(
    Y ~
        c2phonation.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr"),
    data = closure.pl,
    method = "ML"
)

summary(closure.pl.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr")
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.24355    0.14312  22.662   <2e-16 ***
## c2phonation.ord.L  0.05709    0.20260   0.282    0.778    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df       F  p-value    
## s(X)                          7.429  8.356 128.934  < 2e-16 ***
## s(X):c2phonation.ordvoiceless 1.894  2.367   7.644 0.000309 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.523   Deviance explained = 52.7%
## -ML = 4461.8  Scale est. = 29.302    n = 1432
closure.pl.gam.null <- bam(
    Y ~
        s(X, bs = "cr"),
    data = closure.pl,
    method = "ML"
)

compareML(closure.pl.gam, closure.pl.gam.null)
## closure.pl.gam: Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord, 
##     bs = "cr")
## 
## closure.pl.gam.null: Y ~ s(X, bs = "cr")
## 
## Chi-square test of ML scores
## -----
##                 Model    Score Edf Chisq    Df p.value Sig.
## 1 closure.pl.gam.null 4469.756   3                         
## 2      closure.pl.gam 4461.755   6 8.000 3.000   0.001  ** 
## 
## AIC difference: -11.18, model closure.pl.gam has lower AIC.
plot_smooth(closure.pl.gam, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * X : numeric predictor; with 30 values ranging from -14.986800 to 71.049000.

plot_diff(closure.pl.gam, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
##  * X : numeric predictor; with 100 values ranging from -14.986800 to 71.049000.

## 
## X window(s) of significant difference(s):
##  -7.165364 - 11.084655
##  40.632303 - 71.049000
closure.pl.gamm <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
#        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = closure.pl,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(closure.pl.gamm)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.12889    0.40861   7.657 3.91e-14 ***
## c2phonation.ord.L -0.00129    0.54255  -0.002   0.9981    
## c2place.ord.L     -1.03057    0.57581  -1.790   0.0737 .  
## vowel.ord.L       -0.72943    0.55364  -1.318   0.1879    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                   edf  Ref.df      F p-value    
## s(X)                            8.785   8.915 47.967 < 2e-16 ***
## s(X):c2phonation.ordvoiceless   5.262   6.307  1.741 0.21884    
## s(X,rec.date)                 207.831 227.000 88.103 < 2e-16 ***
## s(X):c2place.ordvelar           8.647   8.899 34.922 < 2e-16 ***
## s(X):vowel.ordo                 6.395   7.429  3.651 0.00219 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.987   Deviance explained = 98.9%
## fREML = 2424.9  Scale est. = 0.82697   n = 1432
closure.pl.gamm.null <- bam(
    Y ~
#        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
#        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = closure.pl,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(closure.pl.gamm, closure.pl.gamm.null)
## closure.pl.gamm: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord)
## 
## closure.pl.gamm.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, rec.date, 
##     bs = "fs", xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord)
## 
## Chi-square test of fREML scores
## -----
##                  Model    Score Edf Chisq    Df p.value Sig.
## 1 closure.pl.gamm.null 2430.334  12                         
## 2      closure.pl.gamm 2424.913  15 5.421 3.000   0.013  *  
## 
## AIC difference: -10.11, model closure.pl.gamm has lower AIC.
AIC(closure.pl.gamm, closure.pl.gamm.null)
##                            df      AIC
## closure.pl.gamm      243.6748 4015.338
## closure.pl.gamm.null 240.2176 4025.451
plot_smooth(closure.pl.gamm, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -14.986800 to 71.049000. 
##  * rec.date : factor; set to the value(s): 24/03/2017 14:49:08.

plot_diff(closure.pl.gamm, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -14.986800 to 71.049000. 
##  * rec.date : factor; set to the value(s): 24/03/2017 14:49:08.

## 
## Difference is not significant.
acf_plot(resid(closure.pl.gamm), split_by=list(closure.pl$rec.date))

3.6.1 AR1 model

r1 <- start_value_rho(closure.pl.gamm)

closure.pl.gam.AR <- bam(
    Y ~
        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
#        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = closure.pl,
    method = "fREML",
    rho = r1,
    AR.start = closure.pl$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(closure.pl.gam.AR)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         3.1253     0.4017   7.779 1.57e-14 ***
## c2phonation.ord.L  -0.0335     0.5372  -0.062   0.9503    
## c2place.ord.L      -1.0248     0.5664  -1.809   0.0706 .  
## vowel.ord.L        -0.7378     0.5462  -1.351   0.1770    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                   edf  Ref.df      F p-value    
## s(X)                            8.751   8.908 42.801 < 2e-16 ***
## s(X):c2phonation.ordvoiceless   5.019   6.109  1.417 0.26593    
## s(X,rec.date)                 201.924 227.000 41.057 < 2e-16 ***
## s(X):c2place.ordvelar           8.600   8.896 31.216 < 2e-16 ***
## s(X):vowel.ordo                 6.057   7.195  3.311 0.00231 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.986   Deviance explained = 98.8%
## fREML =   2135  Scale est. = 0.7632    n = 1432
closure.pl.gam.AR.null <- bam(
    Y ~
#        c2phonation.ord +
        c2place.ord +
        vowel.ord +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
#        s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord), bs = "cr",
    data = closure.pl,
    method = "fREML",
    rho = r1,
    AR.start = closure.pl$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(closure.pl.gam.AR, closure.pl.gam.AR.null)
## closure.pl.gam.AR: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") + 
##     s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord)
## 
## closure.pl.gam.AR.null: Y ~ c2place.ord + vowel.ord + s(X, bs = "cr") + s(X, rec.date, 
##     bs = "fs", xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, 
##     bs = "cr") + s(X, by = vowel.ord)
## 
## Chi-square test of fREML scores
## -----
##                    Model    Score Edf Chisq    Df p.value Sig.
## 1 closure.pl.gam.AR.null 2139.752  12                         
## 2      closure.pl.gam.AR 2134.956  15 4.796 3.000   0.022  *
## Warning in compareML(closure.pl.gam.AR, closure.pl.gam.AR.null): AIC is
## not reliable, because an AR1 model is included (rho1 = 0.446313, rho2 =
## 0.446313).
## Warning in compareML(closure.pl.gam.AR, closure.pl.gam.AR.null): Only small difference in fREML...
acf_plot(resid(closure.pl.gam.AR), split_by=list(closure.pl$rec.date))

acf_resid(closure.pl.gam.AR, split_pred = "AR.start")

plot_smooth(closure.pl.gam.AR, view = "X",
            plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 30 values ranging from -14.986800 to 71.049000. 
##  * rec.date : factor; set to the value(s): 24/03/2017 14:49:08.

plot_diff(closure.pl.gam.AR, view = "X",
          comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * X : numeric predictor; with 100 values ranging from -14.986800 to 71.049000. 
##  * rec.date : factor; set to the value(s): 24/03/2017 14:49:08.

## 
## Difference is not significant.